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This paper describes the time-filtered Navier-Stokes (TFNS) approach capable of 
capturing unsteady flow structures important for turbulent mixing in the combustion 
chamber and two different subgrid models used to emulate the major processes oeeurring in 
the turbulence-chemistry interaction. These two subgrid models are termed as LEM-like 
model and EUPDF-like model, respectively. Two-phase turbulent combustion in a single- 
element lean-direct-injeetion (LDI) combustor is calculated by employing the TFNS/LEM- 
like approach as well as the TFNS/EUPDF-like approach. Results obtained from the TFNS 
approach employing these two different subgrid models are compared with each other, along 
with the experimental data, followed by more detailed comparison between the results of an 
updated calculation using the TFNS/LEM-like model and the experimental data. 


Acronym 


CFD 

= computational fluid dynamics 

DNS 

= direct numerical simulation 

LES 

= large-eddy simulation 

RANS 

= Reynolds-averaged Navier-Stokes approach 

URANS 

= unsteady RANS 

TFNS 

= time-filtered Navier-Stokes approach 

SGS 

= subgrid scale 

SFC 

= subfilter component 

FCP 

= filtering-control parameter 

LEM 

= linear-eddy mixing 

EUPDF 

= Eulerian probability density function 

LDI 

= lean-direct injection 

VBB 

= vortex-breakdown bubble 

PVC 

= precession vortex core 

1.0 Introduction 


Together with rig testing and diagnostics, combustion CFD is now a major tool for combustion technology 
development. It is being used to complement, and sometimes to substitute the rig test. This leads to reduced costs, 
deepened insight, and improved foresight. Nevertheless, high-fidelity simulation of liquid combustion in practical 
engineering devices remains an elusive target in spite of significant advances in combustion modeling and 
simulation over the past decade. The main difficulty lies in the fact that combustion operates at the intersection of 
fluid dynamics, fuel chemistry, and multi-phase physics. Consequently, in addition to high-fidelity models for 
turbulence, chemical kinetics, and liquid fuel atomization and transport, models for accurately accounting for their 
interactions also are required. Furthermore, these intricate physical and chemical phenomena are present throughout 
a broad range of time scales and length scales in engineering devices. 

Some of these modeling and simulation issues are discussed in the next section, followed by a description of the 
TFNS formulation for two-phase flow and the subfilter closure models in Section 3. In Section 4, we describe the 
estimation of the filtered, chemical-reaction source terms via two different subgrid models of turbulence-chemistry 
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interaction in the framework of TFNS. In Section 5, we summarize the experimental setup and the computational set 
up employed to assess the TFNS/LEM-like and TFNS/EUPDF-like approaches. This is followed by presenting some 
of the computational results in Section 6 and concluding remarks in Section 7. 

2.0 Liquid-Fueled Turbulent Combustion Modeling Challenges 

The most accurate and straightforward numerical approach to fluid flow problems is to solve the Navier-Stokes 
equations without filtering and approximation, other than numerical discretizations, whose errors are to be 
controlled by highly accurate numerical schemes. This approach is known as the direct numerical simulation (DNS). 
Presently, due to its demand on computing resources, DNS is not yet practical for engineering applications. Besides, 
even when the governing equations are solved directly in DNS, the use of models to accommodate the multiphase 
formulation and interaction is still unavoidable. Modeling and simulation of liquid-fuel injection and spray transport 
is a very difficult task. In fact, many existing large-eddy simulation (LES) of spray flow and combustion invoke the 
same liquid-phase models as those used in the traditional Reynolds-averaged Navier-Stokes (RANS) approach, with 
a subgrid-scale (SGS) model used for the gas-phase turbulence. LES uses spatial filtering to explicitly account for 
flow structures larger than the filter width, the filtering in LES leads to unknown terms in the filtered equations, and 
these so-called SGS terms need to be modeled to form a closed set of the fluid flow governing equations. For 
chemically reacting flow, evaluation of the filtered term arising from the production and destruction of chemical 
species is quite challenging, whether the approach is the traditional RANS or LES, as the turbulence-chemistry 
interaction is present throughout the turbulence spectrum. Furthermore, in the case of liquid-fueled combustion, the 
flame structure is usually very complex and locally (both in space and time) can range from nonpremixed to 
premixed burning. 

In any event, a prerequisite for accurate simulation of turbulent combustion in the combustor is the ability of the 
turbulence model to capture the unsteady turbulent structures responsible for the fuel-oxidizer-products mixing. 
Approaches such as the large-eddy simulation (LES) and the time filtered Navier-Stokes simulation (TFNS) are 
capable of capturing the dynamically important, unsteady turbulent flow structures. In the framework of the 
conventional LES, the filtered equations are established by applying a spatial filter to the exact form of the 
governing equations. The filter width is typically the local grid size; in addition, the eddy viscosity uses the local 
grid size as a model parameter. Therefore, the grid resolution and the model fidelity are formally linked, and, in 
principle, a grid independent solution cannot be reached. An explicitly filtered LES approach can be used to mitigate 
this issue (Ref. 1), and this is an area which needs further development. Last but not the least, the need to use fine 
grids when performing LES is of paramount importance, because, in LES, the subfilter field coincides with the 
subgrid field. 

In the framework of TFNS, the filtered equations are established by applying a temporal filter to the exact form of 
the governing equations. The filter width does not relate to the time step of the numerical solution, and the eddy 
viscosity contains a so called filtering-control parameter, which is nominally defined as the ratio of a (conceptual) 
temporal filter width to a characteristic integral time scale of the turbulent flow. Since the grid resolution and the 
model fidelity are not formally linked, in principle, a grid independent solution can be attained. It should be pointed 
out that TFNS is not LES, nor hybrid RANS/LES, nor, in general, unsteady RANS (URANS). When performing 
TFNS, grids should numerically support the spatial gradients of the filtered variables. Since the smoothness of these 
spatial gradients correlates to the specified value of the filtering-control parameter, numerical accuracy demands that 
the employed grid sizes and grid distribution be consistent with the specified value of the filtering-control parameter. 
This is a reflection of the characteristic of TFNS, i.e., the subfilter field is not the subgrid field. And unlike the scale 
adaptive schemes typically adopted in the hybrid RANS/LES, we do not adapt the eddy viscosity according to the 
local grid resolution. The TFNS approach is summarily presented in Section 3. 

An integral part of turbulent combustion is the turbulence-chemistry interaction. As mentioned above, this 
interaction occurs throughout the entire turbulence spectrum, and the flames usually exhibit multiregime behavior. 
In any approach based on filtered transport equations, whethere it is LES or TFNS or RANS, evaluation of the 
filtered terms arising from the production and destruction of chemical species is quite challenging. To this end, 
various algorithms and models have been reported in the open literature. Here, in Section 4, we present two different 
subgrid models developed for the TFNS approach. 
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Surrogate of the targeted fuel, e.g. the Jet- A, and its transport properties as well as its chemical kinetics should 
also be determined. Fuel chemistry models ranging from global reduced mechanism involving a few species and 
reaction steps to detailed mechanism involving several hundreds of species and reaction steps exist. Global reactions 
are computationally less intensive, but they lack emissions and extinction information. Kinetics calculations using 
detailed mechanism consume long CPU time as well as large computer memory; hence, they are rarely, if ever, 
performed on the flight during the simulation. Typically, libraries of pre-computed tables are established by using 
detailed mechanism, then, linked to the simulation through interpolation procedure. Tabulation strategy such as the 
artificial neural network (ANN) has been used to further speed up the kinetics calculations while reducing the 
memory requirement. 

For liquid-fueled combustion, models for primary atomization, secondary droplet breakup, droplet vaporization, 
and droplet transportation are needed. The performances of these models directly impact the local fuel vapor 
distribution in the computed flow field, and can play a major role in the overall accuracy of the simulation. 

For design and analysis of practical engineering devices, modeling the geometry is as important as modeling the 
physical and chemical processes. Accurate simulation of turbulent combustion in combustors requires precise 
representation of various geometrical features and high-quality mesh distribution in critical regions to numerically 
capture the controlling flow structures and dynamical processes. 

3.0 Mathematical Formulation 

The conservation equations for compressible reacting flow are solved using the TFNS approach. To simulate spray 
combustion, Lagrangian droplet model is concurrently solved with the Eulerian gas phase. In the following, we will 
briefly describe the definitions of time-filtered quantity, the gas-phase equations, and the coupling between the 
gaseous field and the spray field. We will then make general comments on the spray modeling without presenting 
the liquid-phase equations, followed by comments on the geometry modeling. 

3.1 Time-Filtered Quantities </>(x,t) and (f) 

In the case of compressible turbulent flow, we often use two distinct but closely related time-filtered quantities. 
One is denoted by (j) ( X, t ) , which is defined as 


— /• +00 

(j)(x ,0 = f </>(x,t')G(t-t')dt' 

J —00 


( 1 ) 


The integration is over the entire time domain — oo < t ' < + co . Where (f) is an unfiltered turbulent variable and 
G(t-t') is a temporal filter with a constant filter width A r . Furthermore, this temporal filter satisfies the 
following conditions: 
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The other filtered quantity is denoted by (f)(x,t) , which is defined as 

Hx,t) = ^zr 
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(4) 
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That is, (j) (j c,t) is determined by the familiar Favre filtering, which is density-weighted. 


3.2 Gas-Phase TFNS Equations 

Applying a temporal filter with a constant filter width to the exact equations of conservation of mass, momentum, 
energy, and chemical species, as well as the equation of state, the following filtered governing equations in the 
TFN S framework are obtained, 
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where /.t is the filtered mass density, £/. is the filtered velocity vector, i 3 is the filtered pressure determined from 
the filtered equation of state, E is the filtered total energy per unit mass, zv is the filtered viscous stress, q t is the 
filtered heat flux vector, and S, is the source term due to the subfilter kinetic energy. In addition, Y is the filtered 
mass fraction of the m-th species, g . is the filtered mass flux vector of the m-th species. S is the filtered 
reaction-source term of the m-th species, MW is the molecular weight of the m-th species, R is the universal 
gas constant, and T is the filtered temperature, N s is the total number of species. Subscript ‘ / iq ' denotes source 
term due to the liquid phase. Super script ‘ sfc ’ denotes the subfilter component (SFC) which requires closure. 

The filtered viscous stress tensor, the filtered heat flux vector, and the filtered mass flux vector of the m-th species 
are approximated by 
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( 12 ) 
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where Sy is the filtered strain-rate tensor, and h m is the filtered specific enthalpy of the m-th species. Here, the 

molecular viscosity ( // ) is approximated by the Sutherland’s law based on the filtered temperature ( T ). In 
addition, the thermal conductivity ( K ) and the diffusion coefficient of the m-th species are approximated by 

K = ]uc p IP r 


D * 


m 


p Scj 


where C is the specific heat at constant pressure for the gaseous mixture, P r is the Prandtl number, and Sc T is the 
turbulent Schmidt number. 

The SFC terms that require closure are 
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In the present study, we neglect terms such as qf , Ilf , and Tf. The closure for the rest of the SFC terms is 
described below. 

3.2.1 Momentum Transport Closure 

A general constitutive relationship between subfilter turbulent stresses and filtered strain rate (Ref. 2) leads to a 
nonlinear model forz^ , i.e. 
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The model coefficients C , A 3 and A s are constrained by the realizability condition and the rapid distortion 
theory limit. They are formulated as (Ref. 3): 
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The coefficient / is a function of the filtering-control parameter (FCP) which is nominally defined as the ratio of 
a (conceptual) time filter width A r to an integral time scale J 1 , i.e., FCP = A r / T . Furthermore (Ref. 4), 


7 l T ) ( T J i T 


By definition, the value of the parameter FCP and the value of the coefficient f are always between 0 and 1 . 
This model uses the concept of subfilter eddy viscosity which is defined as 

v t = f -Cu-k 2 1 £ 
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Flere, the subfilter turbulent kinetic energy and its dissipation rate ( k, £ ) are supplied by the following model 
equations: 
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where <J k and <J £ are model constants, and the source terms are given by 
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where C £l and C.. 2 arc the model coefficients. We have adopted the commonly used values of C £l = 1.45 and 
C £ 2 = 1.92 in the present work. 

In practice, the content of the subfilter turbulent kinetic energy is regulated by the (specified) filtering-control 
parameter FCP. Since the subfilter turbulent kinetic energy is, by definition, a fraction of the total turbulent energy, 
its value should be smaller than the RANS turbulent kinetic energy. 


3.2.2 Energy ; Transport Closure 

The subfilter total energy flux, Ef c , is also modeled using the subfilter eddy viscosity and gradient assumption, a 
nonlinear model is formulated as (Ref. 3) 
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here, Pr 7 is the turbulent Prandtl number, C u is the specific heat at constant volume for the gaseous mixture. The 
coefficients C, and C 2 are yet to be extensively calibrated. In the current simulations they are set to be C x = C 2 = - 
0.24. The subfilter thermal conductivity ( K T ) is evaluated via 


K r = ptjC / Pr r 


It is noted here that the filtered total energy ( E ) is evaluated as 


E = d + -U i U i 
2 


where e is the filtered specific internal energy of the mixture. The subfilter viscous work is modeled as 
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3.2.3 Species Transport Closure 


Similar to the subfilter total energy flux, the subfilter mass flux of the m-th species is modeled as 


F£ =~pDm T 
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where Sc T is the turbulent Schmidt number. 

In the present study, the filtered reaction rate of the m-th species (i.e. S m in Eq. (8)) is estimated via subgrid 
models of turbulent mixing and combustion, and these models are described in Section 4. 

3.3 Coupling between Gaseous Field and Spray Field 

The coupling between the gaseous and spray field is through the interphase exchange terms. The effects of the 
spray field on the gaseous field are accounted for via the ‘liquid’ source terms in the gas-phase equations (Eq. (5), 
Eq. (6), Eq. (7), and Eq. (8)), and they are modeled as (Ref. 5) 
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where V c is the volume of the computational cell, fl k is the number of droplets in the k-th drop group, th k is the 
vaporization rate of droplets in the k-th group, ll kj is the i-th velocity component of the k-th group, p k is the 
density of droplets in the k-th group, r k is the droplet radius in the k-th group, h t is the specific enthalpy at the 
droplet surface, l k e(/ is the effective latent heat of vaporization of droplets in the k-th group, and £ m is the 
fractional vaporization rate for the m-th species. 



3.4 Liquid-Phase Modeling 


The governing equations of the liquid phase are based on a Lagrangian formulation where the spray particle 
position and velocity are described by a set of ordinary differential equations. Various submodels, such as the 
droplet drag model and the drop vaporization model, are needed to simulate the transport of a vaporizing spray 
particle. The specification of the fuel injection condition plays a major role in the fidelity of the simulation. 
Common practice is to specify the starting droplet condition using correlations of droplet sizes calibrated by relevant 
experimental data. In addition to the use of correlation, various models for primary atomization and secondary 
droplet breakup also have been employed. A more detailed description of the liquid phase modeling and the two- 
way coupling between the liquid-phase and the gas-phase transport can be found in Ref. 5. 

3.5 Geometry Modeling 

To improve the representation of complex geometrical features and to enable high-quality mesh distribution in 
critical regions, we have established capabilities of generating and using arbitrary polyhedral mesh derived from 
hanging-node elements as well as capabilities of adaptive mesh refinement using information from the flow 
solutions. Description and demonstration of these capabilities can be found in Ref. 6. 


4.0 Subgrid Models for Mixing and Combustion 

Simulation of turbulent combustion requires modeling the effects of turbulence-chemistry interaction processes 
occurring throughout the entire turbulence spectrum, and various approaches have been invoked in many previous 
studies. Examples of these approaches include the eddy-break-up (EBU) model (Ref. 7), the thickened flame model 
(Ref. 8), the flamelet-based method (Ref. 9), the conditional moment closure (CMC) method (Ref. 10), the filtered 
mass density function/probability density function (FDF/PDF) method (Ref. 11), and the linear-eddy mixing (LEM) 
model (Ref. 12). Sometimes, a rudimentary model, termed here as the un-mixed (umx) model, is employed. The un- 
mixed model does not involve itself with the mixing processes happening within the CFD computational cells. 
Rather, the (highly nonlinear) filtered reaction source terms are evaluated simply by using the filterd temperature 
and filtered speciese mass fractions directly supplied by the solutions of the filtered governing equations. 

Invoking the FDF/PDF method or the LEM model usually leads to a hybrid algorithm. The velocity and the 
density over the entire computational domain are obtained by solving the filtered Navier-Stokes equations via finite- 
difference scheme or finite-volume scheme. The scalars (i.e., the energy and the species mass fractions) over the 
entire computational domain are obtained by solving the FDF/PDF transport equation or the LEM transport equation 
via stochastical elements such as the Monte Carlo particles (in the case of FDF/PDF method) or the linear cells (in 
the case of LEM model), i.e., the filtered conservation equations for the scalars are not employed. Examples of 
hybrid LES/FDF and hybrid LES/LEM can be found in Ref. Hand Ref. 12, respectively. By their nature, both the 
FDF/PDF method and the LEM model can handle the multiregime flame structures, but the FDF/PDF method 
usually lumps the molecular diffusion and the small scale turbulence effects into a single process known as micro- 
mixing. An Eulerian Monte Carlo implementation of the FDF/PDF method, termed as the EUPDF model, is 
described in relative detail in Ref. 13. Examples of hybrid RANS/EUPDF and hybrid TFNS/EUPDF can be found in 
Ref. 14 and Ref. 15, respectively. 

Under the current TFNS approach, the energy and the species mass fractions over the entire computational domain 
are obtained by solving the filtered conservation equations. The LEM model or the EUPDF model is used as a 
subgrid model for mixing and combustion to provide an estimation of the filtered reaction source terms appearing in 
the filtered conservation equations, and they are termed as the LEM-like model and the EUPDF-like model, 
respectively. More specifically, (i) instead of using the LEM model or the EUPDF model to directly establish the 
entire species field without actually solving the filtered conservation equations of species, the LEM model or the 
EUPDF model is used locally (both in space and time) to only provide approximations of filtered reaction source 
terms in the filtered equations which are now being solved, (ii) instead of solving the unfiltered quantities using the 
LEM equations or the EUPDF equations, the same mathematical forms of the LEM equations or the EUPDF 
equations for the unfiltered quantities are analogously adopted as a subgrid model for mixing and combustion of 
time-filtered scalars, with the molecular transport coefficients substituted by the ‘effective’ transport coefficients. 
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4.1 LEM-like Model 


We have adapted the LEM model described in Ref. 12 into the current TFNS framework, and, for the reasons 
mentioned above, we call it LEM-like. The LEM-like model is implemented in terms of a fractional splitting 
technique; it is divided into a supergrid process and a subgrid process. The supergrid process emulates the 
convection of the scalar field by grid-resolved velocity field across the surfaces of the computational cell. The 
subgrid process, which occurs within each computational cell, consists of four operators: (a) (effective) molecular 
diffusion, (b) finite-rate chemical reaction, (c) volumetric expansion caused by the heat release, and (d) stochastic 
stirring due to the subgrid eddies. 

4.1.1 Supergrid Process 

The macro mixing of the scalars is implemented by a Lagrangian transfer of LEM elements across the surfaces of 
the CFD computational cells. This Lagrangian transport is known as the ‘splicing’, and each LEM element carries its 
own volume, filtered temperature, and filtered species density. Referring to Fig. 1, for example, the outward mass 
through the right side of a cell computed from filtered velocity and density is equivalent to 1 .5 LEM elements and 
colored in red. The outward mass through the bottom side of the mesh is equivalent to 2.5 LEM elements and 
colored in magenta. Similarly, the inward mass through the top and left sides of the mesh are equivalent to 6 LEM 
elements. Splicing will result in 14 LEM elements in this computational cell. In general, splicing will cause different 
computational cell to have different number of LEM elements. To avoid programming complexities in a parallel 
environment, the LEM domain is re-gridded to have the same fixed number of elements, and each element is of the 
same volume. Conservation of mass is maintained during the re-gridding procedure. 

Arrows — Convection direction 
Initial number of lumps - 12 
Light Green - inward mass (1.5) 

Dark Green - inward mass (4.5) 

Red - outward mass (1 .5) 

Magenta - outward mass (2.5) 

Current number of lumps -- 14 

Final number of lumps will become 12 after re-gridding 
via equalizing the volume. 

Figure 1. Schematic illustrating the splicing algorithm used for scalar convection in the LEM-like model. 

4.1.2 Subgrid Process 

Within each CFD computational cell, a one-dimensional domain consisting of a fixed number of LEM elements is 
employed, the governing equations have the following form: 
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where /L stirring and F- stirring represent the subgrid turbulent mixing, and they are accounted for by employing a 
stochastic rearrangement of the LEM elements, known as the triplet mapping, more details about this mapping 
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procedure can be found in Ref. 12. Here, D m „■ is the effective mass diffusion coefficient, /V„ is the effective heat 

conductivity, and d) m (Y l ,Y 2 , - ■ • ,Y N ) indicates that the reaction kinetics of unfiltered species mass fractions is 
formally used for the reaction kinetics of filtered species mass fractions. 

The solution of the above one-dimensional LEM-like equations over NLEN elements gives a subgrid field of Y , 

namely, Y^~ l , Y^~ 2 , • • • , Y l ~ NLEN . Their average, Y 2EM (t + At) , provides a prediction of Y m {t + At) at the center 
of the computational cell. This predicted value is then used to construct an estimaion of the filtered source term 
( S m ) needed to solve Eq. (8). In the current study, the estimation is via 

m At 

where Y^ FNS (t) is the TFNS solution of the (filtered) species mass fraction at the computational cell center. 

4.2 EUPDF-like Model 

We have adapted the EUPDF model described in Ref. 13 into the current TFNS framework, and, for the reasons 
mentioned above, we call it EUPDF-like. The EUPDF-like model is implemented by making use of an operator 
splitting scheme; it is divided into a supergrid process and a subgrid process. The supergrid process accounts for the 
convection as well as the diffusion of the scalar field by grid-resolved velocity field across the surfaces of the 
computational cell. The subgrid process, which occurs within each computational cell, consists of two operators: (a) 
micro mixing, and (b) finite -rate chemical reaction. 

4. 2. 1 Supergrid Process 

The grid-resolved convection and diffusion of the scalar field is implemented in an Eulerian fashion. Referring to 
Figure 2, unlike the LEM splicing algorithm, the Eulerian convection and diffusion process is achieved through the 
content replacement of particles, as each Monte Carlo particle carries the information of the mass fraction of species, 
the enthalpy, the temperature, but not the volume. For example, the outward mass through the right side of a cell, 
computed from the grid-resolved velocity and density; and normalized by the mass in the cell, is equal to 5% and 
colored in red. The outward mass through the bottom side of the mesh is equal to 10% and colored in magenta. 
Similarly, the inward mass through the top and left sides of the mesh are equivalent to 40% and colored in dark and 
light green. If the total number of the particles is 20 in this cell, then one particle (5% of 20) will be randomly 
selected from its right adjacent cell, and its contents are copied to one randomly selected particle in the current cell. 
The same procedure is applied to other particles in colors. For those remaining particles (45% here) the contents are 
randomly shuffled within the cell. 

Arrows - Convection direction 
Light Green - inward mass (15%=3) 

Dark Green - inward mass (25%=5) 

Red - outward mass (5%=1 ) 

Magenta - outward mass (10%=2) 

Number of Monte Carlo particles - 20 
Contents of randomly selected particles are 
replaced by those of randomly selected particles 
from adjacent cells. 



Figure 2 Schematic illustrating the process of the convection and diffusion of the scalar field in the EUPDF- 
like model 
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4.2.2 Subgrid Process 

Within each CFD computational cell, an ensemble of Monte Carlo particles is employed. The number of particles 
within each CFD computational cell is fixed. These particles are notional and, for each particle, the governing 
equations have the following form: 



-cSl{Y n -< Y m >) + d>„ 
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( 28 ) 
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where C . is an empirical constant , £2 is a measure of the turbulence frequency (see, for example, Ref. 15), and < > 
denotes the ensemble average over the Monte Carlo particles in a given CFD computational cell. In addition, 
(O m (Y l ,Y 2 ,---,Y N ) indicates that the reaction kinetics of unfiltered species mass fractions is formally used for the 
reaction kinetics of filtered species mass fractions. 

The solution of the above EUPDF-like equations over NPAR particles gives a subgrid field of Y , 

namely, Y^' , Y ^~ 2 , • • • , Y^~ NPAR . Their average, Y FUPDF ( t + At) , provides a prediction of Y m ( t + At) at the center 
of the computational cell. This predicted value is then used to construct an estimation of the filtered source term 
( S m ) needed to solve Eq. (8). In the current study, the estimation is via 

_ f^PDF {t + At) _ f T FN S {t ) 

m At 

where Y PFNS (t ) is the TFNS solution of the (filtered) species mass fraction at the computational cell center. 


5.0 Experimental and Computational Setup 

The modeling-and-simulation capabilities outlined in the previous sections are embodied in a NASA in-house 
code known as the National Combustion Code (Refs. 16-22). Its target application is for low-emissions, fuel-flexible 
combustors, including the lean-direct- injection (LDI) combustors. 

A typical single-element LDI combustor is illustrated in Fig. 3. More detailed description of the combustor 
geometry and the test rig can be found in Ref. 23. Each element consists of an air passage with an upstream air 
swirler and a converging-diverging venturi section. The fuel is injected through the center of the swirl er and the fuel 
tip is at the throat of the venturi. The air swirler has six helical, axial vanes with downstream vane angles of 60- 
degree. The air then dumps into a combustion chamber with a square cross-section. Velocity measurements were 
taken with a two-component Laser Doppler Velociometry (LDV) system, temperature measurements were taken 
with thermocouples, emissions data were gathered via an isokinetic probe and gas analyzer, and droplet data were 
collected with a Phase Doppler Particle Analyzer (PDPA). The temperature measurements may suffer from radiation 
and conduction loss and other loss due to the presence of liquid droplets. In regions close to the dump plane, the 
gas-phase velocity measurements may suffer from interferences caused by the presence of small liquid particles. 
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The extent of the computational domain is shown in Fig. 3. At the air inflow boundary, the air flow speed is 20.14 
m/s, the density is 1.19 kg per cubic meter, and the static temperature is 294.28 °K. The operating pressure of the 
combustor is approximately 1 atm, and the measured pressure drop (as a percentage of the air inlet pressure) during 
the experiments was measured at 4 %. At the combustor chamber exit, an outlet boundary condition facilitating the 
convection of pressure disturbances out of the computational domain is applied (Ref. 24). Adiabatic and non-slip 
conditions, along with a generalized wall function (Ref. 25), are imposed on the wall surfaces. 

The liquid fuel is injected at 0.415 g/s, which gives a global equivalence ratio of 0.75. The specification of the 
starting condition for the fuel spray is critical for high-fidelity simulations. In this study, the following droplet size 
distribution is used at the injection point (Ref. 5): 


— = 4.21xl0 6 [— ] 35 e 

n d 32 


16 . 98(— ) 04 
^32 


dd 

d 32 


(30) 


where n is the total number of the droplets and dll is the number of droplets in the size range between d and 
d + dd . This correlation requires specifying a Sauter mean diameter. d J2 , as well as the number of droplet 

groups. In the present simulation, the Sauter mean diameter is set to be 50 microns, the number of injected droplet 
groups is set to be 10. At each injection event, a total of 240 spray particles enter into the computational domain at a 
point located in the immediate neighborhood of the injector tip. The magnitude of the droplet injection velocity is 
specified to be 12 m/s. The hollow spray cone has a cone angle of 90 degree with a 30-degree half-cone angle. 
Subsequently, these starting droplets will undergo evaporation, but without experiencing secondary breakup. 

In this study, the liquid fuel Ci 2 H 2 3 is adopted as the surrogate for the experimental Jet-A fuel, and a single-step, 
five-species global reduced mechanism (see e.g. Ref. 26) is employed for the chemical reactions. 


Table 1 Single Step (Global) Chemistry Model. 



Reaction 

A 

(mole - cm - sec- K) 

n 

E 

(cal/mole) 

1 

4 C12H23 + 71 02 => 48 C02 + 46 H20 

GLO/C12H23 0.10/ 

GLO/02 1.65/ 

8.60E+11 

0.00 

3.00E+4 


The computational domain consists of approximately 862,000 hexahedral elements. Current TFNS (with FCP = 
0.255) solutions are obtained by an iteratively implicit, finite-volume scheme that is (nominally) second-order 
accurate in space and time. The LEM-like equations are solved using a standard fractional splitting, finite-difference 
scheme, and 24 LEM elements per CFD computational cell are employed. The EUPDF-like equations are solved by 
an operator splitting scheme, and 24 Monte Carlo particles per CFD computational cell are used. The solution for 
the ordinary differential equations of spray particle position, size, and velocity are advanced via a second-order 
accurate Runge-Kutta method. 

6.0 Results 

While recognizing that the inherent uncertainties in the current models of liquid fuel atomization/injection, fuel 
evaporation, and fuel chemistry often can be overwhelming, this LDI burner and its measured data has been used to 
assess as well as to guide the development of modeling-and-simulation capabilities for two-phase turbulent 
combustion in LDI combustors. Results from the RANS approach were reported in Ref. 26 and Ref. 27. Results 
from the LES/LEM approach were presented in Ref. 12. Results from the LES/flamelet approach were described in 
Ref. 9. And results from TFNS/EUPDF were shown in Ref. 15. In the following, some of the results obtained by the 
TFNS/EUPDF-like approach and the TFNS/LEM-like approach are presented. 
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First of all, it is noted here that the TFNS results of non-reacting flow and their comparison with the experimental 
data can be found in Ref. 28. They reveal the major flow structures in the LDI combustor. Fig. 4 is a snapshot of the 
unsteady flow field. Embedded in this figure are the instantaneous iso-surface of the zero axial velocity component 
colored by the effective eddy viscosity and six instantaneous stream lines emanating from the upstream of the 
swirler, going through the converging-diverging nozzle, then passing through the combustion chamber. Major flow 
structures in the LDI combustor are visualized via the iso-surface of the zero axial velocity and the iso-surface of a 
relatively low pressure. The iso-surface of the zero axial velocity is also known as the vortex-breakdown bubble 
(VBB). The iso-surface of a sufficiently low pressure captures the precession vortex core (PVC). Fig. 5 is a snapshot 
of the PVC and VBB. The dark blue region is a vortex core, which is formed near the venturi throat and extends into 
the combustor chamber. This spiraling vortex core rotates and breaks, it changes in seemingly random way in space 
and time. The light green surfaces are the iso-surfaces of the zero axial velocity. In addition to the VBB, there are 
some small structures near the dump plane and in the corner region. Fig. 6 shows the instantaneous streamlines 
around the PVC, they start from the upstream of the swirler, and yield a complex, seemingly random pattern. Some 
of them spiral around the dark blue surface indicating that the dark blue region is indeed a vortex core. As 
demonstrated in Ref. 12, the dynamics of the PVC and the VBB, as well as their interactions, are critical to the fuel- 
air mixing and the flame stability in the LDI combustors. 

Subsequently, the same grid is used for the two-phase reacting flow calculations, and Fig. 7 is a snapshot of the 
fuel-droplet distribution. The performance of the turbulence-chemistry interaction models is at first highlighted in 
terms of the contours of the computed temperature field, followed by comparisons of averaged temperature and 
averaged velocity at several locations. In addition, time history of flow variables collected at a ‘probe’ location is 
used to illustrate the unsteady features from using different subgrid models. 

The center-plane (z=0) temperature field from the LEM-like model, the EUPDF-like model, and the un-mixed 
(umx) model is compared in Fig. 8. In the case of the un-mixed model and the EUPDF-like model, the filtered 
temperature shown is an average over 10000 time-steps (i.e., 0.01 seconds). In the case of the LEM-like model, the 
filtered temperature shown is an average over 5000 time-steps (i.e., 0.005 seconds). The features of the averaged 
temperature field obtained by using the LEM-like model are closer to what is typically observed in the experiments, 
as compared to the un-mixed and the EUPDF-like models. 

Comparisons of averaged temperature and averaged velocity at several locations are shown from Fig. 9 to Fig. 13. 
They are averages over 50,000 time-steps (i.e., 0.05 seconds). The centerline averaged temperature and the 
centerline averaged axial velocity are given in Fig. 9 and Fig. 10, respectively. Fig. 11 depicts the radial profile of 
the averaged temperature at 10 mm downstream of the dump plane. The radial profiles of the averaged axial 
velocity and the averaged azimuthal velocity at 15 mm downstream of the dump plane are presented in Fig. 12 and 
Fig. 13, respectively. 

Fig. 14, Fig. 15 and Fig. 16 illustrate the time history of the gage pressure, the temperature and the axial velocity 
sampled over 50,000 time-steps at a centerline location which is 7.8 mm downstream of the dump plane, obtained by 
using the EUPDF-like model. Their counterparts obtained by using the LEM-like model are shown in Fig. 17, Fig. 
18 and Fig. 19. It is quite clear that the EUPDF-like model is prone to yield ‘noisy’ results. These noises are 
numerical artifacts caused by the nature of the Monte Carlo process. 

Based on the above comparison, subsequent investigation has been focused on the results from the LEM-like 
model, and new sets of calculations using an updated version of the National Combustion Code (NCC) have been 
carried out. This NCC version contains bug fixes and newly implemented spray-wall interaction models. One set of 
computed data is examined here. The sampling duration of this set is 0.1 second (from time-step 430,001 to time- 
step 530,000), and various time averages are constructed by averaging over these 100,000 time-steps. 

In the following, comparisons between time averages of filtered quantities and experimental mean values are 
presented. The discrepancies are due to a variety of reasons. In addition to the obvious computational reasons such 
as the grid resolution and model fidelity, the measurement errors aggravated by the presence of liquid droplets also 
play an appreciable role. 

The averaged temperature along the center line downstream of the dump plane (located at x=0.0072 m) is shown 
in Fig. 20. Near the dump plane, the significant causes of the discrepancy are the spray injection model, droplet 
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vaporization model, turbulent combustion model, as well as the radiation and conduction loss of the thermocouples. 
Further downstream, the subgrid model for turbulent mixing plays a major role. Fig. 21 presents the averaged axial 
velocity along the combustor center line. In the immediate neighborhood of the dump plane, there was experimental 
difficulty in sorting the seeder particles used in the LDV from the small yet high axial-momentum spray particles. In 
addition, computational heat-release tied to the spray transport and combustion model is another cause of the 
discrepancy. These factors also contribute to the discrepancies between the computational data and the experimental 
data shown in Fig. 22, Fig. 23 and Fig. 24. Fig. 22 is the radial profile of the averaged temperature at 5 mm 
downstream of the dump plane. The computational temperature is the pure gas temperature, however, the measured 
temperature consists of contributions from both pure gas and droplets. Since the sprays distribute as a hollow cone 
donut ring, the discrepancy in temperature is much larger in the donut ring, as compared to the discrepancy in the 
central region. Fig. 23 is the radial profile of the averaged axial velocity at 5 mm downstream of the dump plane. 
The measured gas axial velocity is all positive in the central region in spite of the recirculation zone existing at this 
location. This can be attributed to the presence of the small droplets which have positive velocities but are difficult 
to be separated from the seeding particles used in the LDV system. The radial profile of the averaged azimuthal 
velocity at 5 mm downstream of the dump plane is presented in Fig. 24. The presence of the small droplets which 
have low swirling velocities reduces the measured gas phase swirl velocity in the central region. 

Radial profile of the averaged temperature at 20 mm downstream of the dump plane is shown in Fig. 25. Radial 
profiles of the averaged axial velocity and the averaged azimuthal velocity at 29 mm downstream of the dump plane 
are shown in Fig. 26 and Fig. 27, respectively. There are still discrepancies, it is apparent that grid resolution in the 
cross-section needs to be improved. Radial profiles of averaged temperature, averaged axial velocity, and averaged 
azimuthal velocity at further downstream locations are presented from Fig. 28 to Fig. 33. 

Fig. 34 is the unsteady pressure recorded at a centerline point which is 7.8 mm downstream of the dump plane. 
Fig. 35 and Fig. 36 illustrate the intermittent feature of the mass fraction of the oxygen and the mass fraction of the 
fuel vapor at this ‘probe’ location. In Fig. 37 and Fig. 38, contour plots of time-averaged quantities in the center 
plane (z=0) are presented. Time-average of filtered fuel (C 12 H 23 ) vapor mass fraction is shown in Fig. 37, along with 
two snapshot solutions. Similarly, time-average of filtered temperature, along with two snapshot solutions, are 
shown in Fig. 38. The time-averaged field gives the impression of an orderly, symmetric burning, while the 
snapshots suggest that, in parts of the flame region, the supply of fuel vapor is intermittent, and pockets of 
significant temperature variation exist. 


7.0 Concluding Remarks 

Two-phase turbulent combustion in a single-element LDI combustor has been simulated by employing the 
TFNS/LEM-like approach and the TFNS/EUPDF-like approach. Together with an assumed starting condition of the 
fuel droplets and a single-step, 5-species global reduced kinectics for the Jet-A combustion, computational results 
are obtained by using a computational grid consisting of approximately 862000 hexahedral cells. Results from the 
TFNS/LEM-like and the TFNS/EUPDF-like approaches are compared with each other, together with the 
experimental data. Subsequently, we have focused on the LEM-like model, as it explicitly includes the effects of 
molecular diffusion. These results are useful in gaining deepened insight into the unsteady physical processes in the 
LDI combustor. Nevertheless, in addition to the obvious numerical reasons and other model-related causes, their 
quantitative accuracy often suffers from the inherent uncertainties in the current spray-related models. Future plan 
includes replacing the present single-step combustion chemistry with a multi-step global reduced mechanism, 
achieving a grid-independent solution, and employing atomization as well as secondary breakup models to 
characterize the fuel injection processes. In regions around the injector tip, we also need to find a way to take 
account of the influences of the injected/dense liquid on the flow field of the co-flowing gaseous stream. 
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Geometry of the Single Element 


Grid Distribution for the LDI Combustor (861823 hexahedron elements) 


Figure 3. Swirler geometry and computational domain of a single-element LDI combustor. 



Figure 4. A snapshot of the unsteady flow field. 
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Figure 5. A snapshot of the precessing-vortex core (PVC) and the vortex-breakdown bubble (VBB). 



Figure 6. Instantaneous streamlines around the precession vortex core (PVC). 
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Figure 7. A snapshot of the fuel-droplet distribution in the center plane and a cross-section. 
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Figure 8. Comparison of subgrid models: temperature field in the center plane (z=0). 
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Figure 11. Radial profile of averaged temperature (10 mm downstream of the dump plane). 
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Figure 12. Radial profile of averaged axial velocity (15 mm downstream of the dump plane). 
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Figure 13. Radial profile of averaged azimuthal velocity (15 mm downstream of the dump plane). 
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Figure 14. Pressure history at a center point which is 7.8 mm downstream of the dump plane. 
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Figure 15. Temperature history at a center point which is 7.8 mm downstream of the dump plane. 
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Figure 16. Axial velocity history at a center point which is 7.8 mm downstream of the dump plane. 
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Figure 17. Pressure history at a center point which is 7.8 mm downstream of the dump plane. 



Figure 18. Temperature history at a center point which is 7.8 mm downstream of the dump plane. 
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Figure 19. Axial velocity history at a center point which is 7.8 mm downstream of the dump plane. 



Figure 20. Averaged temperature along the center line. 
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Figure 23. Radial profile of averaged axial velocity (5 mm downstream of the dump plane). 



Figure 24. Radial profile of averaged azimuthal velocity (5 mm downstream of the dump plane). 
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Figure 25. Radial profile of averaged temperature (20 mm downstream of the dump plane). 



Figure 26. Radial profile of averaged axial velocity (29 mm downstream of the dump plane). 












T(k) 



Figure 32. Radial profile of averaged axial velocity (76 mm downstream of the dump plane). 
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Figure 33. Radial profile of averaged azimuthal velocity (76 mm downstream of the dump plane). 



Figure 34. Pressure history at a center point which is 7.8 mm downstream of the dump plane. 
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Figure 35. Oxidizer history at a center point which is 7.8 mm downstream of the dump plane. 



Figure 36. Fuel vapor history at a center point which is 7.8 mm downstream of the dump plane. 
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Figure 37. Fuel vapor distribution in the center plane: time-averaged field and two snapshots. 
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Figure 38. Temperature distribution in the center plane: time-averaged field and two snapshots. 
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